
//！输入nt,nr 得到z=0时的 平面上的数据保存到二维数组中
// double型数据
#ifndef GAUSSIAN_H_
#define GAUSSIAN_H_

#include "./Constants.h"

using std::cout;
using std::endl;

// GaussianDual spatial/temporal pulse focused with a parabolic mirror
class GaussianDual {
private:
    // waist is spatial 1/e^2 radius, duration is full width half maximum (FWHM)
    // focus is positive for converging beams
    double m_z0; 
    double sigma1, sigma2, m_duration, m_Tmax, m_phase0, tau;
    double m_waist1, m_focus1, m_E_peak1, m_wavelength1, m_omega1, m_phase1;
    double m_waist2, m_focus2, m_E_peak2, m_wavelength2, m_omega2, m_phase2;
    double zr1, k01, w_z1, R_z1, eta_z1;//光斑半径//波前曲率半径//Gouy 位移
    double zr2, k02, w_z2, R_z2, eta_z2;//光斑半径//波前曲率半径//Gouy 位移
    
    double radius(double z, double waist, double focus,double zr);
    double curvature(double z, double focus, double zr);
    double gouy(double z, double focus, double zr);

public:
    explicit GaussianDual(double z0, double duration, double T_max, double phase0,
                    double waist1, double focus1, double E_peak1, double wavelength1, double phase1,
                    double waist2, double focus2, double E_peak2, double wavelength2, double phase2)
     :m_z0(z0), m_duration(duration), m_Tmax(T_max), m_phase0(phase0),
     m_waist1(waist1), m_focus1(focus1), m_E_peak1(E_peak1), m_wavelength1(wavelength1),m_phase1(phase1), 
     m_waist2(waist2), m_focus2(focus2), m_E_peak2(E_peak2), m_wavelength2(wavelength2),m_phase2(phase2){ 

        tau = m_Tmax / 2.;
        sigma1 = m_duration / (2. * sqrt(2. * log(2.)));
        sigma2 = sigma1 / 1.;

        zr1 = PI * m_waist1 * m_waist1 / m_wavelength1; //瑞利范围
        k01 = 2. * PI / m_wavelength1;
        m_omega1 = 2. * PI * Constants::c / m_wavelength1;

        w_z1   = radius(m_z0, m_waist1, m_focus1, zr1);
        R_z1   = curvature(m_z0, m_focus1, zr1);
        eta_z1 = gouy(m_z0, m_focus1, zr1);


        zr2 = PI * m_waist2 * m_waist2 / m_wavelength2;
        k02 = 2. * PI / m_wavelength2;
        m_omega2 = 2. * PI * Constants::c / m_wavelength2;
        
        w_z2   = radius(m_z0, m_waist2, m_focus2, zr2);
        R_z2   = curvature(m_z0, m_focus2, zr2);
        eta_z2 = gouy(m_z0, m_focus2, zr2);   
    }

    virtual ~GaussianDual() {

    }
 
    double Electric1(double r,double t);
    double Electric2(double r,double t);
    void operator()(int row, int column, double *r_array, double *t_array, double **E_array);

    double Intensity(double z, double r);//光强空间分布


};

double GaussianDual::radius(double z, double waist, double focus, double zr){
    return waist * std::sqrt(1.0 + std::pow((z - focus) / zr,2));
}

double GaussianDual::curvature(double z, double focus, double zr){
    return (z - focus) * (1.0 + std::pow(zr / (z - focus),2));
}

double GaussianDual::gouy(double z, double focus, double zr){
    return std::atan2(z - focus, zr);
}

double GaussianDual::Electric1(double r, double t) {
        double arg = 0;
        arg += m_omega1 * t + m_phase1;
        arg += - k01 * std::pow(r,2) / (2. * R_z1) ;
        arg += eta_z1 ;
        
        double Er_envelop =  (m_waist1 / w_z1) * std::exp(- pow(r / w_z1, 2)) ;
        double Et_envelop =  std::exp(- 0.5* pow((t - tau) / sigma1, 2));
        double EE = m_E_peak1 * Er_envelop  * Et_envelop * std::cos(arg); 
        return EE;
}

double GaussianDual::Electric2(double r, double t) {
        double arg = 0;
        arg += m_omega2 * t + m_phase2 + m_phase0;
        arg += - k02 * std::pow(r,2) / (2. * R_z2) ;
        arg += eta_z2 ;
        
        double Er_envelop =  (m_waist2 / w_z2) * std::exp(- pow(r / w_z2, 2)) ;
        double Et_envelop =  std::exp(- 0.5* pow((t - tau) / sigma2, 2));
        double EE = m_E_peak2 * Er_envelop  * Et_envelop * std::cos(arg); 
        return EE;
}


void GaussianDual::operator()(int row, int column, double *r_array, double *t_array, double **E_array){
    for(int num_row = 0; num_row < row; num_row++){
        double ri  = r_array[num_row];
        for(int num_column = 0; num_column < column; num_column++){
            double tj   = t_array[num_column];
            E_array[num_row][num_column] = Electric1(ri, tj) + Electric2(ri,tj);
            
        }
    }
}
//输出强度

double GaussianDual::Intensity(double z, double r){
    double wz1 = radius(z, m_waist1, m_focus1, zr1);
    double Erz_envelop1 = (m_waist1 / wz1) * std::exp(-pow(r/wz1, 2));

    double wz2 = radius(z, m_waist2, m_focus2, zr2);
    double Erz_envelop2 = (m_waist2 / wz2) * std::exp(-pow(r/wz2, 2));
    double Amplitude = m_E_peak1 * Erz_envelop1 + m_E_peak2 * Erz_envelop2;

    return Amplitude;
}

// double GaussianDual::Intensity(double z, double r, double t){
//     double arg = 0.0;
//     arg += m_omega1 * t;
//     arg += - k01 * std::pow(r, 2) / (2 * curvature(z));
//     arg += gouy(z);

//     double wz1 = radius(z);
//     double Er_envelop = (m_waist1 / wz1) * std::exp(-pow(r/wz1, 2)) ;

//     double Et_envelop = std::exp(- 0.5* pow((t - tau) / sigma1, 2));

//     double Amplitude = m_E_peak1 * Er_envelop  * Et_envelop * cos(arg);

//     return Amplitude;
// }

//!-----------------------------------------------------------------------------------------------------
// GaussianSingle spatial/temporal pulse focused with a parabolic mirror
class GaussianSingle {
private:
    // 实验上：聚焦后，所有频率的曲率半径、光束半径、gouy相位都等于中心频率的值
    // waist is spatial 1/e^2 radius, duration is full width half maximum (FWHM)
    // focus is positive for converging beams
    double m_z0, m_duration, m_Tmax, m_phase0;
    double m_waist, m_focus, m_E_peak, m_wavelength, m_omega;
    double zr, k0, tau ;
    double sigma;
    double w_z   ;   //光斑半径
    double R_z   ;   //波前曲率半径
    double eta_z ;   //Gouy 位移  

    double radius(double z);
    double curvature(double z);
    double gouy(double z);

public:
    explicit GaussianSingle(double z0, double duration, double T_max, double phase0,
                            double waist, double focus, double E_peak, double wavelength)
     :m_z0(z0), m_duration(duration), m_Tmax(T_max), m_phase0(phase0), 
     m_waist(waist), m_focus(focus), m_E_peak(E_peak), m_wavelength(wavelength){ 
        
        zr = PI * m_waist * m_waist / m_wavelength; //瑞利范围
        k0 = 2. * PI / m_wavelength;
        m_omega = 2. * PI * Constants::c / m_wavelength;
        tau = m_Tmax / 2.;
        sigma = m_duration / (2. * sqrt(2. * log(2.)));

        w_z = radius(m_z0);
        R_z = curvature(m_z0);
        eta_z = gouy(m_z0);   
    }

    virtual ~GaussianSingle() {

    }
 
    double Electric(double r,double t);
    void operator()(int row, int column, double *r_array, double *t_array, double **E_array);

    double Intensity(double z, double r);//光强空间分布


};

double GaussianSingle::radius(double z){
    return m_waist * std::sqrt(1.0 + std::pow((z - m_focus) / zr,2));
}

double GaussianSingle::curvature(double z){
    return (z - m_focus) * (1.0 + std::pow(zr / (z - m_focus),2));
}

double GaussianSingle::gouy(double z){
    return std::atan2(z - m_focus, zr);
}

double GaussianSingle::Electric(double r, double t) {
        double arg = 0;
        arg += m_omega * t + m_phase0;
        arg += - k0 * std::pow(r,2) / (2 * R_z) ;
        arg += eta_z ;
        
        double Er_envelop =  (m_waist / w_z) * std::exp(- pow(r / w_z, 2)) ;

        double Et_envelop =  std::exp(- 0.5* pow((t - tau) / sigma, 2));
        
        double EE = m_E_peak * Er_envelop  * Et_envelop * std::cos(arg); 
     
        return EE;
}


void GaussianSingle::operator()(int row, int column, double *r_array, double *t_array, double **E_array){
    for(int num_row = 0; num_row < row; num_row++){
        double ri  = r_array[num_row];
        for(int num_column = 0; num_column < column; num_column++){
            double tj   = t_array[num_column];
            E_array[num_row][num_column] += Electric(ri, tj);
        }
    }
}
//输出强度

double GaussianSingle::Intensity(double z, double r){
    double wz = radius(z);
    double Erz_envelop = (m_waist / wz) * std::exp(-pow(r/wz, 2));
    double Amplitude = m_E_peak * Erz_envelop ;

    return Amplitude;
}


// double GaussianSingle::Intensity(double z, double r, double t){
//     double arg = 0.0;
//     arg += m_omega * t;
//     arg += - k0 * std::pow(r, 2) / (2 * curvature(z));
//     arg += gouy(z);

//     double wz = radius(z);
//     double Er_envelop = (m_waist / wz) * std::exp(-pow(r/wz, 2)) ;

//     double Et_envelop = std::exp(- 0.5* pow((t - tau) / sigma, 2));

//     double Amplitude = m_E_peak * Er_envelop  * Et_envelop * cos(arg);

//     return Amplitude;
// }

#endif //GAUSSIAN_H_

//!laser propagation(github)(未看懂)
// double df = m_focus1 / (1 + std::pow(m_focus1/zr, 2)); 
// double GaussianSingle::radius(double z)const{
//     double b = std::pow(z/zr,2);
//     if (m_focus1 == 0.0){
//         return m_waist1 * std::sqrt(1+b);
//     }
//     else{
//         double a = std::pow(1-z/m_focus1,2);
//         return m_waist1 * std::sqrt(a+b);
//     }
// }

// double GaussianSingle::curvature(double z)const{
//     if(m_focus1 ==0){
//         return z / (std::pow(z, 2) + std::pow(zr, 2));
//     }
//     else{
//         double a = z - df;
//         return a / (std::pow(a, 2) + df * (m_focus1 - df));
//     }
// }

// double GaussianSingle::gouy(double z)const{
//     if (m_focus1 == 0){
//         return std::atan2(z,zr);
//     }
//     else{
//         double a = z - df;
//         double b = std::sqrt(m_focus1 * df - std::pow(df, 2));
//         return std::atan2(a, b);
//     }
// }